{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import numpy as np  # 矩阵操作\n",
    "import pandas as pd # SQL数据处理\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt   #画图\n",
    "from sklearn.metrics import r2_score  #评价回归预测模型的性能\n",
    "from sklearn.model_selection import train_test_split # 数据分割\n",
    "# 数据标准化\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "# 线性回归\n",
    "from sklearn.linear_model import LinearRegression\n",
    "# 线性模型，随机梯度下降优化模型参数\n",
    "from sklearn.linear_model import SGDRegressor\n",
    "#岭回归 --> L2正则\n",
    "from sklearn.linear_model import  RidgeCV\n",
    "#### Lasso --> L1正则\n",
    "from sklearn.linear_model import LassoCV\n",
    "\n",
    "# 图形出现在Notebook里而不是新窗口\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 读取数据文件并总体概览\n",
    "data_path = '/home/ai/study/week1/exercises/basic/'\n",
    "data_filename = 'day.csv'\n",
    "data_full_filename = data_path + data_filename\n",
    "\n",
    "if not os.path.exists(data_full_filename):\n",
    "    print '[-] file(%s) is not found!', data_full_filename\n",
    "    sys.exit(-1)\n",
    "    \n",
    "data = pd.read_csv(data_full_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 根据数据探索中剔除一些数据\n",
    "# 作业不要求y列 ‘casual'和'registered'，先剔除\n",
    "# 日期列'dteday'非数字，并且'holiday' 'weekday' 'workingday'已经足以表示该列，也剔除\n",
    "# 执行下面代码，重新执行整个notebook，后读入数据已经没有这两列了.....\n",
    "try:\n",
    "  #data = data.drop('dteday', axis = 1)\n",
    "  #data = data.drop('casual', axis = 1)\n",
    "  #data = data.drop('registered', axis = 1)\n",
    "  data = data.drop(['dteday', 'casual', 'registered', 'instant'], axis = 1)\n",
    "except:\n",
    "  pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 准备数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_train = data[data.yr == 0] # 2011年数据用于训练\n",
    "data_test = data[data.yr == 1]  # 2012年数据用于测试\n",
    "\n",
    "# 风速达到高点部分有一部分骑行，先剔除这几个离群点\n",
    "data_train = data_train[data_train.windspeed < 0.5]\n",
    "\n",
    "y_train = data_train['cnt'].values\n",
    "X_train = data_train.drop('cnt', axis = 1)\n",
    "\n",
    "y_test = data_test['cnt'].values\n",
    "X_test = data_test.drop('cnt', axis = 1)\n",
    "\n",
    "columns = X_test.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(364, 11)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(366, 11)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 数据标准化\n",
    "ss_X = StandardScaler()\n",
    "ss_y = StandardScaler()\n",
    "\n",
    "X_train = ss_X.fit_transform(X_train)\n",
    "X_test  = ss_X.transform(X_test)\n",
    "\n",
    "y_train = ss_y.fit_transform(y_train.reshape(-1, 1))\n",
    "y_test  = ss_y.transform(y_test.reshape(-1, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 线性回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>columns</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>[0.629841694384535]</td>\n",
       "      <td>temp</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[0.2950038141022198]</td>\n",
       "      <td>season</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[0.03690725066076217]</td>\n",
       "      <td>weekday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>[0.0065142549413874234]</td>\n",
       "      <td>atemp</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>[0.004933089112578709]</td>\n",
       "      <td>workingday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[0.003712357146231554]</td>\n",
       "      <td>mnth</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[-8.326672684688674e-16]</td>\n",
       "      <td>yr</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[-0.03669761849370721]</td>\n",
       "      <td>holiday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>[-0.06276419978882453]</td>\n",
       "      <td>hum</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>[-0.11594761649624319]</td>\n",
       "      <td>windspeed</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>[-0.22491174014925078]</td>\n",
       "      <td>weathersit</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        coef     columns\n",
       "7        [0.629841694384535]        temp\n",
       "0       [0.2950038141022198]      season\n",
       "4      [0.03690725066076217]     weekday\n",
       "8    [0.0065142549413874234]       atemp\n",
       "5     [0.004933089112578709]  workingday\n",
       "2     [0.003712357146231554]        mnth\n",
       "1   [-8.326672684688674e-16]          yr\n",
       "3     [-0.03669761849370721]     holiday\n",
       "9     [-0.06276419978882453]         hum\n",
       "10    [-0.11594761649624319]   windspeed\n",
       "6     [-0.22491174014925078]  weathersit"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 线性回归模型训练和预测和评估\n",
    "lr = LinearRegression()\n",
    "lr.fit(X_train, y_train)\n",
    "y_test_pred_lr = lr.predict(X_test)\n",
    "y_train_pred_lr = lr.predict(X_train)\n",
    "fs = pd.DataFrame({\"columns\":list(columns), \"coef\":list((lr.coef_.T))})\n",
    "fs.sort_values(by=['coef'],ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The r2 score of LinearRegression on test is -0.6925610091966108\n",
      "The r2 score of LinearRegression on train is 0.7580954504877622\n"
     ]
    }
   ],
   "source": [
    "#测试集\n",
    "print 'The r2 score of LinearRegression on test is', r2_score(y_test, y_test_pred_lr)\n",
    "#训练集\n",
    "print 'The r2 score of LinearRegression on train is', r2_score(y_train, y_train_pred_lr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# L2正则 --> 岭回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The r2 score of RidgeCV on test is -0.6971286301843629\n",
      "The r2 score of RidgeCV on train is 0.7573638050835745\n"
     ]
    }
   ],
   "source": [
    "# 岭回归模型训练和预测和评估\n",
    "\n",
    "#设置超参数（正则参数）范围\n",
    "alphas_v = [ 0.01, 0.1, 1, 10,100]\n",
    "ridge = RidgeCV(alphas = alphas_v, store_cv_values = True)\n",
    "\n",
    "ridge.fit(X_train, y_train)\n",
    "y_test_pred_ridge = ridge.predict(X_test)\n",
    "y_train_pred_ridge = ridge.predict(X_train)\n",
    "print 'The r2 score of RidgeCV on test is', r2_score(y_test, y_test_pred_ridge)\n",
    "print 'The r2 score of RidgeCV on train is', r2_score(y_train, y_train_pred_ridge)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt8VfWZ7/HPQy7c74RbQrjITSyXQEBtvSK2eCkXFbFTe2prx94809ZjO8yLHttxTk9bOdM5Z3psq7adqWc6FYKAqFBFRG07tbJDuAiIIgp7J0Aicg8QSJ7zx16he9KdCyQreyf5vl8vXqy91m+t/awl5snvt9Z6fubuiIiIXKxOqQ5ARETaNiUSERFpFiUSERFpFiUSERFpFiUSERFpFiUSERFpFiUSERFpFiUSERFpFiUSERFplsxUB9AaBgwY4CNGjEh1GCIibUpxcfEH7p7TWLsOkUhGjBhBJBJJdRgiIm2Kme1tSjsNbYmISLMokYiISLOEmkjMbLaZ7TKz3Wa2KMn2B8xsh5ltNbP1ZjY8YVu+mb1oZjuDNiOC9TeY2SYz22xmvzez0WGeg4iINCy0RGJmGcCjwE3ABOBTZjahTrMSoNDdJwHLgUcStj0JLHH3S4EZQHmw/qfAp919CvDvwLfDOgcREWlcmD2SGcBud9/j7lXAU8DcxAbuvsHdK4OPrwN5AEHCyXT3dUG7EwntHOgVLPcGykI8BxERaUSYT23lAtGEzzHg8gba3wusDZbHAkfMbAUwEngJWOTu1cAXgDVmdgo4BlyR7GBmdh9wH0B+fn4zTkNERBoSZo/EkqxLOh2jmd0NFAJLglWZwNXAg8B0YBRwT7DtG8DN7p4H/Avwo2THdPfH3b3Q3Qtzchp9DFpERC5SmIkkBgxL+JxHkmEoM5sFLAbmuPuZhH1LgmGxc8AqYKqZ5QCT3f1PQbulwEfDOgERkbYqdriSH/72LcqPnw79u8JMJBuBMWY20syygbuA1YkNzKwAeIx4Eimvs2/fIHEAzAR2AIeB3mY2Nlh/I7AzxHMQEWmTiiIxfvbqu5ytTjoQ1KJCu0fi7ufM7H7gBSAD+KW7bzezh4GIu68mPpTVAygyM4B97j7H3avN7EFgvcU3FANPBMf8a+BpM6shnlg+H9Y5iIi0RTU1zvLiGFeNHkBun66hf1+oJVLcfQ2wps66hxKWZzWw7zpgUpL1K4GVLRimiEi78od3P6D0yCkW3TS+Vb5Pb7aLiLQzyyIx+nTL4uOXDWqV71MiERFpR45UVvHC9gPMm5JL58yMVvlOJRIRkXbkmc1lVJ2r4c7CYY03biFKJCIi7cjSjVE+ktuLCUN7Nd64hSiRiIi0E2+WHmXH/mMsbMXeCCiRiIi0G8siUbIzOzFncm6rfq8SiYhIO3D6bDWrSkq56SOD6d0tq1W/W4lERKQdeGH7AY6dPteqN9lrKZGIiLQDyyJR8vp25cpR/Vv9u5VIRETauOiHlfxh9yEWTBtGp07JCq+HS4lERKSNKyqOYQZ3FOal5PuVSERE2rDqGmd5JMrVY3JapUBjMkokIiJt2B92f0DZ0dPcmaLeCCiRiIi0acsiUfp0y+LGCa1ToDEZJRIRkTbq8MkqXtx+sFULNCajRCIi0kY9s7mUqurWLdCYTKiJxMxmm9kuM9ttZouSbH/AzHaY2VYzW29mwxO25ZvZi2a2M2gzIlj/OzPbHPwpM7NVYZ6DiEg6cneWRmJMzO3dqgUakwktkZhZBvAocBMwAfiUmU2o06wEKHT3ScBy4JGEbU8CS9z9UmAGUA7g7le7+xR3nwL8EVgR1jmIiKSrN0uPsXP/Me6cntreCITbI5kB7Hb3Pe5eBTwFzE1s4O4b3L0y+Pg6kAcQJJzMYLpd3P1EQjuCNj2BmYB6JCLS4SyLROmc2Yk5k4emOpRQE0kuEE34HAvW1edeYG2wPBY4YmYrzKzEzJYEPZxE84H17n4s2cHM7D4zi5hZpKKi4iJPQUQk/Zw+W82qzUGBxq6tW6AxmTATSbL39D1pQ7O7gUJgSbAqE7gaeBCYDowC7qmz26eA39T35e7+uLsXunthTk7OhUUuIpLGXth+gOMpKtCYTJiJJAYknmUeUFa3kZnNAhYDc9z9TMK+JcGw2Dniw1dTE/bpT3zo7PmQYhcRSVtLN0YZ1q8rV6SgQGMyYSaSjcAYMxtpZtnAXcDqxAZmVgA8RjyJlNfZt6+Z1XYlZgI7ErYvAJ5z99OhRS8ikoaiH1byH++mrkBjMqElkqAncT/wArATWObu283sYTObEzRbAvQAioLHeVcH+1YTH9Zab2bbiA+TPZFw+LtoYFhLRKS9KopE4wUap6WuJEpdmWEe3N3XAGvqrHsoYXlWA/uuAybVs+26FgpRRKTNqK5xlhfHuGZMDkNTVKAxGb3ZLiLSRvz+fIHG9LjJXkuJRESkjVgWidK3WxazJgxMdSj/iRKJiEgbcPhkFeu2H2ReQWoLNCajRCIi0gasCgo0LkyDkih1KZGIiKQ5d2fpxiiT8nozfnBqCzQmo0QiIpLmtpUe5a0Dx9PuJnstJRIRkTRXW6Dxk2lQoDEZJRIRkTR2+mw1z2wu4+aJQ9KiQGMySiQiImnst2/GCzQuKEyfN9nrUiIREUljSzdGye/XjStGpkeBxmSUSERE0tS+Q5X8cc8h7izMS5sCjckokYiIpKmi4iidDG5PowKNySiRiIikofMFGsfmMKR3+hRoTEaJREQkDf3unQr2p2GBxmSUSERE0lBRJEa/7tnMunRQqkNplBKJiEia+fBkFS/uOMC8KblkZ6b/j+lQIzSz2Wa2y8x2m9miJNsfMLMdZrbVzNab2fCEbflm9qKZ7QzajAjWm5l9z8zeDrb9TZjnICLS2laWlHK22tOyQGMyoc2QaGYZwKPAjUAM2Ghmq909ce71EqDQ3SvN7MvAI8DCYNuTwPfcfZ2Z9QBqgvX3AMOA8e5eY2bpVZhfRKQZ3J2iSJTJeb0ZN7hnqsNpkjB7JDOA3e6+x92rgKeAuYkN3H2Du1cGH18H8gDMbAKQGUy3i7ufSGj3ZeBhd68JtpWHeA4iIq1qaywo0NhGeiMQbiLJBaIJn2PBuvrcC6wNlscCR8xshZmVmNmSoIcDcAmw0MwiZrbWzMYkO5iZ3Re0iVRUVDTzVEREWseySJQuWelboDGZMBNJstcwPWlDs7uBQmBJsCoTuBp4EJgOjCI+pAXQGTjt7oXAE8Avkx3T3R9390J3L8zJybnYcxARaTWnqqpZvbmMmz8yhF5d0rNAYzJhJpIY8XsZtfKAsrqNzGwWsBiY4+5nEvYtCYbFzgGrgKkJ254OllcCk0KIXUSk1f12+36OnznHgjbw7kiiMBPJRmCMmY00s2zgLmB1YgMzKwAeI55Eyuvs29fMarsSM4Ham/Srgs8A1wJvhxS/iEirWroxyvD+3bhiVL9Uh3JBQkskQU/ifuAFYCewzN23m9nDZjYnaLYE6AEUmdlmM1sd7FtNfFhrvZltIz5M9kSwzw+A24P13we+ENY5iIi0lr2HTvL6ng+5s3AYZulboDGZ0B7/BXD3NcCaOuseSlie1cC+60gybOXuR4BbWjBMEZGUK4rE4gUap6Z3gcZk0v+VSRGRdq62QOO1Y3MY3LtLqsO5YEokIiIp9to7FRw41jYKNCajRCIikmLLNkbp1z2bG9pAgcZklEhERFLo0IkzvLTzIPML2kaBxmTaZtQiIu1EbYHGtjqsBUokIiIp4+4si0SZPKxPmynQmIwSiYhIimyJHeXtgydY2IZ7I6BEIiKSMrUFGm+dPCTVoTSLEomISAqcqqrm2c1l3DyxbRVoTEaJREQkBda+GS/Q2JZvstdSIhERSYGlG6OM6N+Ny0e2rQKNySiRiIi0svc/OMmf3vuQBW2wQGMySiQiIq2sqDjaZgs0JqNEIiLSimoLNF43bmCbLNCYjBKJiEgreu3tCg4eO8Odhe2jNwJKJCIirWrpxij9u2czc3zbLNCYTKiJxMxmm9kuM9ttZouSbH/AzHaY2VYzW29mwxO25ZvZi2a2M2gzIlj/r2b2XjCj4mYzmxLmOYiItJT2UKAxmdDOxMwygEeBm4AJwKfMbEKdZiVAobtPApYDjyRsexJY4u6XAjOAxDndv+nuU4I/m8M6BxGRlrSypJRzNc6d09v+uyOJwkyJM4Dd7r7H3auAp4C5iQ3cfYO7VwYfXwfyAIKEkxlMt4u7n0hoJyLS5rg7SzdGmTKsD2MHtd0CjcmEmUhygWjC51iwrj73AmuD5bHAETNbYWYlZrYk6OHU+l4wHPZPZtY52cHM7D4zi5hZpKKiojnnISLSbJujR3in/AQL21lvBMJNJMnesvGkDc3uBgqBJcGqTOBq4EFgOjAKuCfY9nfA+GB9P+Bvkx3T3R9390J3L8zJybnIUxARaRnLIjG6ZmVw66S2XaAxmTATSQxITL15QFndRmY2C1gMzHH3Mwn7lgTDYueAVcBUAHff73FngH8hPoQmIpK2KqvO8eyWeIHGnm28QGMyYSaSjcAYMxtpZtnAXcDqxAZmVgA8RjyJlNfZt6+Z1XYlZgI7gn2GBH8bMA94M8RzEBFptrXbDnDizLl2OawF8SGkULj7OTO7H3gByAB+6e7bzexhIOLuq4kPZfUAioJ6M/vcfY67V5vZg8D6IGEUA08Eh/51kGAM2Ax8KaxzEBFpCUsjUUYO6M70EX1THUooQkskAO6+BlhTZ91DCcuzGth3HTApyfqZLRmjiEiY3vvgJG+89yHfmj2uXRRoTKb9vBEjIpKGiiLtq0BjMkokIiIhOVddw/LiGNePG8igXu2jQGMySiQiIiF57Z0Kyo+fYUE7mAWxIUokIiIhWboxyoAe2dxw6cBUhxIqJRIRkRB8cOIM63eWM78gl6yM9v2jtn2fnYhIiqzcFBRobOfDWqBEIiLS4tydZZEoBfl9GNPOCjQmo0QiItLCSmoLNHaA3ggokYiItLiiSJSuWRnc0g4LNCajRCIi0oLiBRr3c8uk9lmgMRklEhGRFrSmnRdoTEaJRESkBS3bGGXUgO4UDm+fBRqTUSIREWkheypO8Mb7H7KgcFi7LdCYTJMTiZldZWafC5ZzzGxkeGGJiLQ9RcUxMjoZt09taFbx9qdJicTMvkN8Stu/C1ZlAf8WVlAiIm3Nueoani6Ocf24HAa24wKNyTS1RzIfmAOcBHD3MqD9v2UjItJEr77dMQo0JtPURFLl7g44gJl1b8pOZjbbzHaZ2W4zW5Rk+wNmtsPMtprZejMbnrAt38xeNLOdQZsRdfb9sZmdaGL8IiKhqi3QOHN8+y7QmExTE8kyM3sM6GNmfw28xJ+nvk3KzDKAR4GbgAnAp8xsQp1mJUChu08ClgOPJGx7Elji7pcCM4Dzc7qbWSHQp4mxi4iEquL4GV5+q5zbpua1+wKNyTTpjN39fxH/Qf80MA54yN1/3MhuM4Dd7r7H3auAp4C5dY67wd0rg4+vA3kAQcLJDKbbxd1P1LYLEtQS4FtNiV1EJGwrS2JBgcb2OwtiQ5p6s7078LK7f5N4T6SrmTX2ymYuEE34HAvW1edeYG2wPBY4YmYrzKzEzJYECQTgfmC1u+9vJOb7zCxiZpGKiopGQhURuTjxAo0xpub3YfTAjnnruKl9sNeAzmaWS3xY63PAvzayT7KHqD1pQ7O7gULiPQ2ATOBq4EFgOjAKuMfMhgILgMZ6Q7j74+5e6O6FOTk5jTUXEbkom/YdYXf5iQ71JntdTU0kFgwt3Qb82N3nE7/v0ZAYkHhl84Cyvziw2SxgMTDH3c8k7FsSDIudA1YBU4ECYDSw28zeB7qZ2e4mnoOISIsrikTplp3BLZOGpjqUlGlyIjGzK4FPA88H6zIb2WcjMMbMRppZNnAXsLrOQQuAx4gnkfI6+/Y1s9quxExgh7s/7+6D3X2Eu48AKt19dBPPQUSkRZ08c45nt5Rxy8Qh9Ojc2I/E9qupieRrwCJghbtvD95qf7mhHYKexP3AC8BOYFmw78NmNidotgToARSZ2WYzWx3sW018WGu9mW0jPkzW4FNiIiKt7flt+zlZVd2hh7Wg8V5FrUqghvgjvHcT/8Ge9H5HIndfA6yps+6hhOVZDey7DpjUyPF7NBaDiEhYiiJRRuV0Z1oHKtCYTFMTya+J9xDeJJ5QREQ6tHcrTrDx/cMsuml8hyrQmExTE0mFuz8baiQiIm1IUSReoPG2DlagMZmmJpLvmNnPgfVA7ZNVuPuKUKISEUlj56preHpTjOvHDWRgz45VoDGZpiaSzwHjiVf9rR3ackCJREQ6nFd2VVBx/EyHfZO9rqYmksnuPjHUSERE2oilkSgDenTm+g5YoDGZpj7++3qSgosiIh1O+fHTvPxWObdPze2QBRqTaWqP5Crgs2b2HvF7JAZ4ULVXRKTDWLmplOoa75DzjtSnqYlkdqhRiIi0AfECjVGmDe/L6IF6ja1WkxKJu+8NOxARkXS3ad9h3q04ySO3X5LqUNKKBvhERJpo2cZYUKBxSKpDSStKJCIiTXDyzDme21rGrZOG0L0DF2hMRolERKQJnt+qAo31USIREWmCZUGBxqn5HbtAYzJKJCIijdhdfoLI3sMsLBzW4Qs0JqNEIiLSiKLiKBmdjPkq0JiUEomISAPOVtfwdHEpM8erQGN9Qk0kZjbbzHaZ2W4zW5Rk+wNmtsPMtprZejMbnrAt38xeNLOdQZsRwfpfmNmWYJ/lZqa3gkQkNK/squCDE2e4U2+y1yu0RGJmGcCjwE3ABOKzK9at11UCFAalVpYDjyRsexJY4u6XAjOA2jndv+Huk4N99hGfzldEJBRLN0bJ6dmZ68flpDqUtBVmj2QGsNvd97h7FfAUMDexgbtvcPfK4OPrQB5AkHAyg+l2cfcTte3c/VjQxoCuNGHKXxGRi1F+/DQbdpVz29RcMlWgsV5hXplcIJrwORasq8+9wNpgeSxwxMxWmFmJmS0JejgAmNm/AAeIz5Hy42QHM7P7zCxiZpGKiormnIeIdFArggKNGtZqWJiJJNkzckl7D2Z2N1AILAlWZQJXE58nfjowCrjn/EHcPwcMBXYCC5Md090fd/dCdy/MyVGXVEQuTG2BxsLhfbkkR7diGxJmIokBiWk8Dyir28jMZgGLgTnufiZh35JgWOwcsAqYmrifu1cDS4HbQ4hdRDq44r2H2VNxkjv1JnujwkwkG4ExZjbSzLKBu4DViQ3MrAB4jHgSKa+zb18zq+1KzAR2WNzoYF8DPgm8FeI5iEgHtXRjlO7ZGdwyUQUaGxNa5TF3P2dm9wMvABnAL919u5k9DETcfTXxoaweQFHwtug+d5/j7tVm9iCwPkgYxcATxIfLfmVmvYLlLcCXwzoHEemYTpw5x/Pb9vPJSUNVoLEJQr1C7r4GWFNn3UMJy7Ma2HcdkGwGxo+1WIAiIkk8v7WMyqpqDWs1kZ5nExGpY1kkxiU53Zma3yfVobQJSiQiIgl2lx+neO9hFk5XgcamUiIREUlQFImR2cmYX5CX6lDaDCUSEZHA2eoant4UY+b4geT07JzqcNoMJRIRkcCGt8r54ESV3mS/QEokIiKBZZF4gcbrVKDxgiiRiIgA5cdOs2FXBbdPzVOBxgukN20aUBSJUnrkFH27ZdOnW9Zf/N2jc6ae6hBpJ54+X6BRN9kvlBJJA17YfpCXdh6sd3tWhtG7azZ96yaZ7sHfXbPo0y3Y3j2+vU/XbLIz9duOSDpxd4oiUWaM6McoFWi8YEokDfj5Zws5V13DkVNnOVJZxeHKsxypPMvhyqqEz1UcPhlft/dQJZujRzhSeZaq6pp6j9ujc2bSHs75pJOwvjYx9VTvRyQ0kb2H2fPBSb583SWpDqVNUiJpRGZGJwb06MyAHk1/FNDdqayq5sipsxw+WfUXySe+HP/7cOVZ9n1YyZHKsxw9dbb+ODrZ+WSTrKfTN/jcJ0g+tcvq/Yg07nyBxkkq0HgxlEhCYGZ075xJ986Z5Pbp2uT9qmuco6cSks7Jv0w6R4LPscOVvFkaX3/mXP29n27ZGUnv8ZxPOt3rJJ+u2fTskkmnTur9SMdw4sw5nt+6n7lThtItWz8SL4auWhrJ6GT0655Nv+7ZF7TfqarqINHUSTonq+K9ooT1pUdOcbiyiqOnzuL1TFKc0cno3TUrSU/nPyed3onDb92y6JKVkfyAImnsuS1lnDqrAo3NoUTSDnTNzqBrdleGXmDv59ipv+zpHE7yd+mR02wvO8bhyipOn62/99M1K4ORA7ozZ8pQ5k4ZypDeTY9HJFWWRaKMHtiDgmEq0HixlEg6qIxORt/u2fS9wN7P6bPVSZPNkcr4/aBN+w7zg7Vv8cPfvsUVI/szvyCX2RMH06tLVkhnInLxdpcfZ9O+Iyy++VI9zNIMSiRyQbpkZTCkd9cGext7D51kVUkZqzaX8q2nt/LtZ97kxksHMa8gl2vH5ugBAEkby2oLNE7NTXUobVqoicTMZgP/h/gMiT939x/U2f4A8AXgHFABfN7d9wbb8oGfE5/33YGb3f19M/s1UAicBd4Avuju9T/uJK1ueP/ufG3WGP7mhtFsiR1lVUkpz24p4/lt++nTLYtbJg5hfkEu04b31W+BkjJnq2tYsSnGDZcOvKCnMuUvhZZIzCwDeBS4EYgBG81stbvvSGhWAhS6e6WZfRl4BFgYbHsS+J67rzOzHkDt4PyvgbuD5X8nnoh+GtZ5yMUzM6YM68OUYX1YfMul/P6dD1i1uZSnN8X49Z/2MaxfV+ZNyWXulFxGD9RLYNK6XlaBxhYTZo9kBrDb3fcAmNlTwFzgfCJx9w0J7V8nSBBmNgHIDKbbxd1PJOxzfupeM3sDUD2DNiAroxPXjx/I9eMHcuLMOV7cfoCVJaU8umE3P355N5PyejNvSi6fnDxU5bulVSzbGGVgz85cO1YFGpsrzESSC0QTPseAyxtofy+wNlgeCxwxsxXASOAlYJG7V9c2NrMs4DPA15IdzMzuA+4DyM/Pv8hTkDD06JzJbVPzuG1qHuXHTrN6S/x+ysPP7eB7a3Zy1egBzC/I5eOXDdJz/RKKg8dOs2FXOV+89hIVaGwBYf5fmmzwO+mbC2Z2N/H7HtcGqzKBq4ECYB+wFLgH+EXCbj8BXnP33yU7prs/DjwOUFhYWM8bE5JqA3t14QtXj+ILV4/inYPHWbW5lFUlZXx96Wa6ZWfwicsGM68gl49d0l//w0uLeXpTjBpHw1otJMxEEiN+o7xWHlBWt5GZzQIWA9e6+5mEfUsShsVWAVcQJBIz+w6QA3wxtOil1Y0Z1JNvfmI8/+3GcUT2HmZlSSnPby1jZUkpA3p0Zs7kocwrGMrE3N66SS8XLV6gMcaMkf0YOaB7qsNpF8JMJBuBMWY2EigF7gL+KrGBmRUAjwGz3b28zr59zSzH3SuAmUAk2OcLwCeAG9y9/rfjpM3q1MmYMbIfM0b247tzJrDhrQpWlZTyb6/v5Zd/eI9ROd2ZPyWXeQW5DOvXLdXhShuz8f3DvPfBSb56/ehUh9JumNdXJ6MlDm52M/C/iT/++0t3/56ZPQxE3H21mb0ETAT2B7vsc/c5wb43Av9IfIisGLjP3avM7BywFzge7LPC3R9uKI7CwkKPRCItfXrSyo5WnmXNm/tZWVLKG+99CEDh8L7MK8jllolDLvjlSumY/tuyLbyw/QBvLL5B9+AaYWbF7l7YaLswE0m6UCJpf2KHK3lmcxmrSkp5p/wEWRnGdeMGMr8gl5njB6rulyR1/PRZZnxvPfMKhvL92yalOpy019REonQsbVJe32589frRfOW6S9ix/xirSkp5ZnMZ63YcpGfnTG6eOIR5BblcPrKfKhnLec9t3R8v0Kib7C1KiUTaNDPjsqG9uWxobxbddCl/fPcQK0tKeW5rGUsjUYb07sLcKbnML8hl3OCeqQ5XUmxZJMqYgT2YogKNLUqJRNqNjE7GVWMGcNWYAfyPeR9h3c6DrCop5Ynf7eFnr77LpUN6MW/KUOaoMnGH9M7B45TsO8K3b1GBxpamRCLtUtfsDOZMHsqcyUP54MQZnt8av0n//bVv8YPfvsWVo/ozryCXmz4ymJ6qTNwhLItEyexkzCtQgcaWppvt0qG898FJVpWUsmpzKXsPVdI5sxOzJgxi/pRcrlFl4nar6lwNV35/PdNH9ONnn5mW6nDaDN1sF0li5IDufOPGsXx91hhKokd4pqSUZ7fu5/mt++nbLYtbJw1lXkEuU/P7aPijHXn5rXIOnazizukqzRcGJRLpkMyMqfl9mZrfl2/fOoHfvVPBypIylkWi/L/X95LfrxvzCnKZN2Uoo3JUmbitWxaJMqhXZ64ZowKNYVAikQ4vK6MTM8cPYub4QRw/fZYXtsdv0v/45Xf45/XvMDmvN/MLcrl18lDNW9EGHTh6mld2lfMlFWgMjRKJSIKeXbK4Y1oed0zL48DR0zy7JV7r67vP7uAfnt/J1WPilYlvnKDKxG2FCjSGTzfbRZpg14F4ZeJnSkopO3qabtkZzA4qE39UlYnTlrtz/f96hUG9urD0i1emOpw2RzfbRVrQuME9+dvZ4/nmx8fxxvsfsqqklOe37WdFSSk5PeOViecX5HLZ0F66SZ9G3njvQ94/VMl/nTkm1aG0a0okIhegUyfjilH9uWJUf7475zI2vFXOypJSnvzj+/zi9+8xemAP5hfkMmfyUFUmTgNLI1F6BCVzJDxKJCIXqUtWBjdNHMJNE4dwpLKK57ftZ1VJKUte2MWSF3YxY0Q/5hYM5ZaJQ+jTTZWJW9vx02dZs20/8wvy6JqtIp5hUiIRaQF9umXz6cuH8+nLhxP9sJLVW8pYsSnG4pVv8t3V27l+3EBum5rLdeNUmbi1PLtlP6fP1rBwum6yh02JRKSFDev358rE28uOsTKoTPzijoP06vLnysQzRqgycZiWRaKMHdSDyXm9Ux1Ku6dEIhISM+Mjub35SG5v/u6m8fzHu4dYVVLK6i1lPLUxytDeXZhbkMsd0/K4RC89tqi3Dx5nc1QFGltLqM8smtlsM9tlZrvNbFGS7Q+Y2Q4z22pm680nAI/+AAAQ3ElEQVRseMK2fDN70cx2Bm1GBOvvD47nZjYgzPhFWkpmRieuGZvDjxZOIfLtWfyfu6YwZlBPHnv1XW74x1e57Sd/4Ddv7OP46bOpDrVdWLYxSlaGMV8FGltFaO+RmFkG8DZwIxAjPg/7p9x9R0Kb64E/uXulmX0ZuM7dFwbbXgG+5+7rzKwHUBO0KwAOA68Ahe7+QWOx6D0SSVflx06zoqSUokiUdytO0iWrE7MvG8yCwmFcOaq/hr4uQtW5Gq74/nouH9mPn96tAo3NkQ7vkcwAdrv7niCgp4C5wPlE4u4bEtq/DtwdtJ0AZLr7uqDdiYR9SoI2IYYu0joG9urCl669hC9eM4rN0SMUFcd4dksZqzaXkdunK7dPy+OOqXnk99ejxE21fudBPjxZpTfZW1GYiSQXiCZ8jgGXN9D+XmBtsDwWOGJmK4CRwEvAInevDiNQkVQzMwry+1KQ35eHbp3AC9sPsLw4dr7e1+Uj+7GgcBg3fWQw3Tvr1mZDlkWiDO7VhWvGqkBjawnzX2SyLkPScTQzuxsoBK4NVmUCVwMFwD5gKXAP8Ismf7nZfcB9APn5+U3dTSTlumRlMHdKLnOn5FJ25BQrNsVYXhzjwaItfOeZN7l54hAWFA5j+oi+6pnXceDoaV59u4KvXDeaDA0LtpowE0kMSOxb5gFldRuZ2SxgMXCtu59J2LckYVhsFXAFF5BI3P1x4HGI3yO5mBMQSbWhfbpy/8wxfPX60UT2HqYoEuX5rfspKo4xvH837piax+3T8hjaR1MHw58LNC4o1LwjrSnMRLIRGGNmI4FS4C7grxIbBDfOHwNmu3t5nX37mlmOu1cAMwHdLZcOy8yYPqIf00f04zufvIy1bx5geXGUf1z3Nj966W2uGj2AO6bl8YnLBnfYFx5rapxlkShXjOrH8P7dUx1OhxJaInH3c2Z2P/ACkAH80t23m9nDQMTdVwNLgB5AUdBF3+fuc9y92sweBNZbfEMx8ASAmf0N8C1gMLDVzNa4+xfCOg+RdNO9c+b5Uvf7DlWyfFOMp4tjfO2pzfTsksknJw/ljml5FAzrWLM8vvH+h+w9VMnXblCBxtamMvIi7UBNjfP6nkMUFcdY+2a8NMjogT24Y1oetxXkMrBXl1SHGLoHlm5m3Y6DvLF4lmprtZCmPv6rRCLSzhw/ffb8fZTivYfpZHDt2BwWFA7jhksH0jmz/f2QPXb6LDO+9xK3Tc3jf86fmOpw2o10eI9ERFKgZ5cs7pqRz10z8tlTcYLlxTFWbCrlK7/eRJ9uWcybEi/L0p7mTnl2S1m8QKPeHUkJ9UhEOoDqGud371SwvDjGizsOUnWuhvGDe7KgcBjzpgylfxufi37uo3/gdFU1v/361e0mOaYD9UhE5LyMTsZ14wZy3biBHK08y+otpSwvjvEPz+3g+2t2MnP8QBYUDuO6cTlktbFpg3cdOM6W6BH++60TlERSRIlEpIPp3S2Lz1w5gs9cOYJdB46zvDjKypJSXtxxkAE9spk3JZcFhcMYN7hnqkNtkqUq0JhyGtoSEc5W1/DqrgqKiqOs31nOuRpnUl5v7piWx5zJQ9N2hseqczVc/j9f4spL+vOTT6tAY0vT0JaINFlWRidmTRjErAmDOHTiDM9sLqOoOMZDz2znfzy3kxsvG8Qd0/K4ZkxOWpUeeWnnQQ5XnmWBbrKnlBKJiPwn/Xt05vNXjeTzV43kzdKjLC+O8czmUp7fup9BvTpz29S8tJmMa1kkypDeXbhmjAo0ppISiYjU6/wMjzeP5+Wd5SwvjvH4a3v46SvvMjW/DwsKh3HrpCH07JLV6rHtP3qK196u4KvXq0BjqimRiEijOmdmcNPEIdw0cQjlx06zsqSUouIYf7diG3//7PaUTMb1dHFQoHGahrVSTYlERC7IwF5d+OK1l3DfNaPYEjtKUSTK6laejCteoDHGlaP6a9KvNKBEIiIXxcyYMqwPU4b14b/fOoEXdxykKBJtlcm4/vTeh+z7sJJv3KgCjelAj/+KSIsqO3IqPvQVifL+oUq6Z2e0+GRc31i6mZd2HmTj4lkdtmx+a9DjvyKSEkP7dOWr14/mK9ddQvHewxRFYjy3tazFJuM6dvosa7bt545peUoiaUKJRERCYWYUjuhH4Yh+fGfOBNZui89D39zJuFZvLuPMuRoWTtdN9nShRCIioeuWncnt0+I9keiHlSwvjs9DnzgZ14JpeUxpwmRcRZEo4wf3ZGJu71aKXhoTanU2M5ttZrvMbLeZLUqy/QEz22FmW81svZkNT9iWb2YvmtnOoM2IYP1IM/uTmb1jZkvNLD1rN4hIUsP6deMbN47ld9+6nn//68u58dJBrNgUY/5P/oMb/+k1fvbqu5QfO51037cOHGNL7Ch3Fg5TgcY0EloiMbMM4FHgJmAC8Ckzm1CnWQlQ6O6TgOXAIwnbngSWuPulwAygdk73HwL/5O5jgMPAvWGdg4iEp1Mn46OXDOBHC6ewcfEsfnj7RPp0zeIHa9/iiu+v53P/8gZrtu3nzLnq8/vUFmicpwKNaSXMoa0ZwG533wNgZk8Bc4EdtQ3cfUNC+9eBu4O2E4BMd18XtDsRrDdgJvBXwT6/Ar4L/DTE8xCRkPXsksXC6fksnB6fjOvpTTGeLv7zZFxzJw9lXkEuq0pK+fiEwfTrroGIdBJmIskFogmfY8DlDbS/F1gbLI8FjpjZCmAk8BKwCOgLHHH3cwnH1K8mIu3IqJwefPMT43ngxnH8fvcHLC+O8ZuNUX71x70ALCjMS3GEUleYiSTZAGbSl1bM7G6gELg2WJUJXA0UAPuApcA9wOoLOOZ9wH0A+fn5FxC2iKSDjE7GtWNzuHZsTnwyrq1llB4+xdUq0Jh2wkwkMSDx+bw8oKxuIzObBSwGrnX3Mwn7liQMi60CrgB+CfQxs8ygV5L0mADu/jjwOMRfSGyRMxKRlOjdLYvPXDG88YaSEmE+tbURGBM8ZZUN3EWdHoWZFQCPAXPcvbzOvn3NrPZXj5nADo+/hr8BuCNY/1ngmRDPQUREGhFaIgl6DPcDLwA7gWXuvt3MHjazOUGzJUAPoMjMNpvZ6mDfauBBYL2ZbSM+TPZEsM/fAg+Y2W6gP/CLsM5BREQap1pbIiKSVFNrbYX6QqKIiLR/SiQiItIsSiQiItIsSiQiItIsSiQiItIsHeKpLTOrAPZe5O4DgA9aMJyWorgujOK6MIrrwrTXuIa7e6OlBDpEImkOM4s05fG31qa4LoziujCK68J09Lg0tCUiIs2iRCIiIs2iRNK4x1MdQD0U14VRXBdGcV2YDh2X7pGIiEizqEciIiLNokRSh5ktMbO3zGyrma00sz71tJttZrvMbLeZLWqFuBaY2XYzqzGzep/CMLP3zWxbUE059EqVFxBXa1+vfma2zszeCf7uW0+76uBana8+HVI8DZ6/mXU2s6XB9j+Z2YiwYrnAuO4xs4qEa/SFVorrl2ZWbmZv1rPdzOyfg7i3mtnUNIjpOjM7mnCtHgo7puB7h5nZBjPbGfy/+LUkbcK9Xu6uPwl/gI8Tny8e4IfAD5O0yQDeBUYB2cAWYELIcV0KjANeAQobaPc+MKAVr1ejcaXoej0CLAqWFyX77xhsO9EK16jR8we+AvwsWL4LWJomcd0D/N/W+veU8L3XAFOBN+vZfjPxqbmN+KR3f0qDmK4DnkvBtRoCTA2WewJvJ/nvGOr1Uo+kDnd/0f88J/zrxGdhrGsGsNvd97h7FfAUMDfkuHa6+64wv+NiNDGuVr9ewfF/FSz/CpgX8vc1pCnnnxjvcuAGM0s2XXVrx5US7v4a8GEDTeYCT3rc68RnTh2S4phSwt33u/umYPk48fmfcus0C/V6KZE07PPEs3hduUA04XOMv/wPlyoOvGhmxcG89ekgFddrkLvvh/j/aMDAetp1MbOImb1uZmElm6ac//k2wS8yR4lP3Bampv53uT0YDlluZsOSbE+FdP1/8Eoz22Jma83sstb+8mBItAD4U51NoV6vMOdsT1tm9hIwOMmmxe7+TNBmMXAO+HWyQyRZ1+zH35oSVxN8zN3LzGwgsM7M3gp+k0plXK1+vS7gMPnB9RoFvGxm29z93ebGVkdTzj+Ua9SIpnzns8Bv3P2MmX2JeK9pZshxNUUqrldjNhEvKXLCzG4GVgFjWuvLzawH8DTwdXc/Vndzkl1a7Hp1yETi7rMa2m5mnwVuBW7wYICxjhiQ+JtZHlAWdlxNPEZZ8He5ma0kPnzRrETSAnG1+vUys4NmNsTd9wdd+PJ6jlF7vfaY2SvEf5tr6UTSlPOvbRMzs0ygN+EPozQal7sfSvj4BPH7hukglH9TzZH4w9vd15jZT8xsgLuHXoPLzLKIJ5Ffu/uKJE1CvV4a2qrDzGYTnxd+jrtX1tNsIzDGzEaaWTbxm6OhPfHTVGbW3cx61i4Tf3Ag6RMmrSwV12s18Nlg+bPAX/SczKyvmXUOlgcAHwN2hBBLU84/Md47gJfr+SWmVeOqM44+h/j4ezpYDfyX4GmkK4CjtUOZqWJmg2vva5nZDOI/Xw81vFeLfK8BvwB2uvuP6mkW7vVq7ScM0v0PsJv4WOLm4E/tkzRDgTUJ7W4m/nTEu8SHeMKOaz7x3yrOAAeBF+rGRfzpmy3Bn+3pEleKrld/YD3wTvB3v2B9IfDzYPmjwLbgem0D7g0xnr84f+Bh4r+wAHQBioJ/f28Ao8K+Rk2M6/vBv6UtwAZgfCvF9RtgP3A2+Pd1L/Al4EvBdgMeDeLeRgNPMrZiTPcnXKvXgY+20rW6ivgw1daEn1s3t+b10pvtIiLSLBraEhGRZlEiERGRZlEiERGRZlEiERGRZlEiERGRZlEiEWmAmZ1o5v7Lg7fmG2rzijVQObmpbeq0zzGz3za1vUhzKJGIhCSotZTh7nta+7vdvQLYb2Yfa+3vlo5HiUSkCYI3gpeY2ZsWn+9lYbC+U1AKY7uZPWdma8zsjmC3T5PwRr2Z/TQoELndzP6+nu85YWb/aGabzGy9meUkbF5gZm+Y2dtmdnXQfoSZ/S5ov8nMPprQflUQg0iolEhEmuY2YAowGZgFLAnKh9wGjAAmAl8ArkzY52NAccLnxe5eCEwCrjWzSUm+pzuwyd2nAq8C30nYlunuM4CvJ6wvB24M2i8E/jmhfQS4+sJPVeTCdMiijSIX4SriVXCrgYNm9iowPVhf5O41wAEz25CwzxCgIuHznUFp/8xg2wTiZS0S1QBLg+V/AxIL8NUuFxNPXgBZwP81sylANTA2oX058VI1IqFSIhFpmvommWpo8qlTxGtoYWYjgQeB6e5+2Mz+tXZbIxJrGJ0J/q7mz//vfoN4jbPJxEcYTie07xLEIBIqDW2JNM1rwEIzywjuW1xDvLji74lP/NTJzAYRn2611k5gdLDcCzgJHA3a3VTP93QiXv0X4K+C4zekN7A/6BF9hvj0ubXGkh7Vn6WdU49EpGlWEr//sYV4L+Fb7n7AzJ4GbiD+A/tt4jPTHQ32eZ54YnnJ3beYWQnx6rB7gD/U8z0ngcvMrDg4zsJG4voJ8LSZLSBenfdkwrbrgxhEQqXqvyLNZGY9PD4rXn/ivZSPBUmmK/Ef7h8L7q005Vgn3L1HC8X1GjDX3Q+3xPFE6qMeiUjzPWdmfYBs4B/c/QCAu58ys+8Qnxt7X2sGFAy//UhJRFqDeiQiItIsutkuIiLNokQiIiLNokQiIiLNokQiIiLNokQiIiLNokQiIiLN8v8BY9bQVydNfMsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8499da20d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('alpha is:', 10.0)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef_lr</th>\n",
       "      <th>coef_ridge</th>\n",
       "      <th>columns</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>[0.629841694384535]</td>\n",
       "      <td>[0.34826087095216535]</td>\n",
       "      <td>temp</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[0.2950038141022198]</td>\n",
       "      <td>[0.2755046061268229]</td>\n",
       "      <td>season</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[0.03690725066076217]</td>\n",
       "      <td>[0.03562222196325582]</td>\n",
       "      <td>weekday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>[0.0065142549413874234]</td>\n",
       "      <td>[0.28444804233012116]</td>\n",
       "      <td>atemp</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>[0.004933089112578709]</td>\n",
       "      <td>[0.004812629670809032]</td>\n",
       "      <td>workingday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[0.003712357146231554]</td>\n",
       "      <td>[0.018903727196517617]</td>\n",
       "      <td>mnth</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[-8.326672684688674e-16]</td>\n",
       "      <td>[0.0]</td>\n",
       "      <td>yr</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[-0.03669761849370721]</td>\n",
       "      <td>[-0.03452196424442788]</td>\n",
       "      <td>holiday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>[-0.06276419978882453]</td>\n",
       "      <td>[-0.06622928339553169]</td>\n",
       "      <td>hum</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>[-0.11594761649624319]</td>\n",
       "      <td>[-0.10962147088244117]</td>\n",
       "      <td>windspeed</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>[-0.22491174014925078]</td>\n",
       "      <td>[-0.21585168109657585]</td>\n",
       "      <td>weathersit</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     coef_lr              coef_ridge     columns\n",
       "7        [0.629841694384535]   [0.34826087095216535]        temp\n",
       "0       [0.2950038141022198]    [0.2755046061268229]      season\n",
       "4      [0.03690725066076217]   [0.03562222196325582]     weekday\n",
       "8    [0.0065142549413874234]   [0.28444804233012116]       atemp\n",
       "5     [0.004933089112578709]  [0.004812629670809032]  workingday\n",
       "2     [0.003712357146231554]  [0.018903727196517617]        mnth\n",
       "1   [-8.326672684688674e-16]                   [0.0]          yr\n",
       "3     [-0.03669761849370721]  [-0.03452196424442788]     holiday\n",
       "9     [-0.06276419978882453]  [-0.06622928339553169]         hum\n",
       "10    [-0.11594761649624319]  [-0.10962147088244117]   windspeed\n",
       "6     [-0.22491174014925078]  [-0.21585168109657585]  weathersit"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 这里是原样抄写的\n",
    "mse_mean = np.mean(ridge.cv_values_, axis = 0)\n",
    "plt.plot(np.log10(alphas_v), mse_mean.reshape(len(alphas_v),1)) \n",
    "\n",
    "#这是为了标出最佳参数的位置，不是必须\n",
    "#plt.plot(np.log10(ridge.alpha_)*np.ones(3), [0.28, 0.29, 0.30])\n",
    "\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show()\n",
    "\n",
    "print ('alpha is:', ridge.alpha_)\n",
    "\n",
    "# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性\n",
    "fs = pd.DataFrame({\"columns\":list(columns), \"coef_lr\":list((lr.coef_.T)), \"coef_ridge\":list((ridge.coef_.T))})\n",
    "fs.sort_values(by=['coef_lr'],ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# L1正则 --> Lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ai/anaconda2/lib/python2.7/site-packages/sklearn/linear_model/coordinate_descent.py:1094: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The r2 score of LassoCV on test is -0.8121199460602377\n",
      "The r2 score of LassoCV on train is 0.7236939499716459\n"
     ]
    }
   ],
   "source": [
    "# Lasso模型训练和预测和评估\n",
    "\n",
    "#############################################################################\n",
    "# 如果不区分年份，把所有数据测试\n",
    "# 数据不进行归一化：\n",
    "# 如果Lasso不设置超参数，测试集上r2score: 0.4920926666976039\n",
    "# 如果设置如下超参数, 测试集上r2score: 0.8130081405966603\n",
    "# 设置超参赛，lasso-alpha为0.01, ridge-alpha为0.1, 是因为0.1 * 0.1 == 0.01 ?\n",
    "#\n",
    "# 数据进行归一化：\n",
    "# lasso-alpha: 0.01, ridge-alpha: 10.0\n",
    "# lasso-r2score on test: 0.8206380871933154, 相比不归一化有提升\n",
    "#############################################################################\n",
    "\n",
    "alphas_v = [ 0.01, 0.1, 1, 10,100]\n",
    "\n",
    "#lasso = LassoCV()\n",
    "lasso = LassoCV(alphas = alphas_v)\n",
    "\n",
    "lasso.fit(X_train, y_train)\n",
    "y_test_pred_lasso = lasso.predict(X_test)\n",
    "y_train_pred_lasso = lasso.predict(X_train)\n",
    "print 'The r2 score of LassoCV on test is', r2_score(y_test, y_test_pred_lasso)\n",
    "print 'The r2 score of LassoCV on train is', r2_score(y_train, y_train_pred_lasso)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8VfWd//HX52YlLGFJ2AKICwKyQ7DWvWot7qICUsXujJ2202W6zfiY2mnnNzOOnelMax2kHX+O1mJwrdal7lLbWnPDGmRH5SYBErZAAmT9zh/3ECMm4Ybk3HOX9/PxyINLzjf3vHMh951zvveerznnEBERAQgFHUBERBKHSkFERNqoFEREpI1KQURE2qgURESkjUpBRETaqBRERKSNb6VgZvebWbWZlXey/Ttmttr7KDezFjMb7FceERE5MfPrzWtmdiFQBzzonJt8grHXAN90zl3iSxgREYlJpl937JxbYWZjYxy+EFgWy8CCggI3dmysdysiIgBlZWV7nHOFJxrnWynEyszygDnAV2MZP3bsWMLhsL+hRERSjJm9H8u4RJhovgb4o3NuX2cDzGyxmYXNLFxTUxPHaCIi6SURSuFmTnDqyDm31DlX7JwrLiw84dGPiIicpEBLwczygYuA3waZQ0REonybUzCzZcDFQIGZVQB3AlkAzrkl3rC5wIvOuXq/coiISOz8fPXRwhjGPAA84FcGERHpnkSYUxARkQShUhARkTaBv09BRBKTc47mVkdzi6OxpZXmllaaWhxNLa00tbTS3HrstqO5pdUb0+5zra1tt5uO39bSSlOrAy0H3C3FYwdz4Zn+vgJTpSDiI+fcB0+QzY4m74my/RPkh55km6NPls0tH31CPfbE3NzqPvIE3HTsSdu7jw/GdPCk3Damoyfs9k/o/j9hm/m+i5Ry+0WnqxREksHeugZuXvoWe+sbP3iSb4n+pu0nM8jKCJEVMjIzQtHbGUZWRojMDCPb+zMzFCI7I0R2Zoi+HxoT/dpj47M62pYZIjNkZGeGyAx9+H4/ND7U0b7bb4/eV5Y3LiNkmFoh4agURHrBY2UVbKmu4+bZo+mTndH2ZPnhJ8oQ2RkffvJuv/34J+bMUIjsTG9MZvsn/g/GZ4T0pCq9S6Ug0kPOOUrCEYpPGcS/3jg16DgiPaJXH4n0UPj9/WyvqWfB7NFBRxHpMZWCSA+VlEbol5PJVVNHBB1FpMdUCiI9cOhoE8+u3ck100aQl62zsZL8VAoiPfDMmp0caWphwewxQUcR6RUqBZEeKAlHGD+sP9NG5QcdRaRXqBRETtLGXQdZEznA/Nmj9Xp7SRkqBZGTVFIaITsjxNwZRUFHEek1KgWRk9DQ3MKTqyr55KRhDO6bHXQckV6jUhA5CS+u382Bw00sKNZ7EyS1qBRETsLycISigX04/4yCoKOI9CqVgkg3RfYd5s2te5hXPIqQrj0kKUalINJNj5ZVADBPp44kBakURLqhpdXxWDjCBeMKKRrYJ+g4Ir1OpSDSDW9u3UNV7VFNMEvKUimIdENJ6Q4G983msrOGBh1FxBcqBZEY7a1r4KV3djN3RhE5mRlBxxHxhW+lYGb3m1m1mZV3MeZiM1ttZuvN7A2/soj0hidXVdLU4rRugqQ0P48UHgDmdLbRzAYC9wLXOucmAfN8zCLSI845SkojzBgzkDOH9Q86johvfCsF59wKYF8XQz4NPOGc2+GNr/Yri0hPrYocYEt1nSaYJeUFOadwJjDIzF43szIzuy3ALCJdWl4aIS87g6unjQw6ioivglwqKhOYBVwK9AH+bGZvOec2Hz/QzBYDiwHGjNFiJhJf9Q3NPLOmiqunjqBfjlZXk9QW5JFCBfCCc67eObcHWAFM62igc26pc67YOVdcWFgY15Aiz67dSX1jiyaYJS0EWQq/BS4ws0wzywM+BmwIMI9Ih0rCEU4v7MvMMYOCjiLiO9+Ohc1sGXAxUGBmFcCdQBaAc26Jc26Dmb0ArAVagV855zp9+apIELZWH6Ls/f3cceVEra4macG3UnDOLYxhzN3A3X5lEOmpktIImSFj7kytribpQe9oFulEY3Mrj6+s5LKJwyjolxN0HJG4UCmIdOKVDbvZV9/IgrM1wSzpQ6Ug0omScIQR+blcOE6veJP0oVIQ6UDVgSO8sbmGm2aNIkOrq0kaUSmIdOCxsgqcg/m6rIWkGZWCyHFaWx3LwxHOO2MIowfnBR1HJK5UCiLH+dO2vVTsP8KC2bqkiqQflYLIcUrCEfL7ZHH5WcOCjiISdyoFkXb21zfy+/JdzJ1RRG6WVleT9KNSEGnnqdWVNLa06uJ3krZUCiKeY6urTR2Vz8QRA4KOIxIIlYKIZ11lLRt3HdLLUCWtqRREPI+URsjNCnHtdK2uJulLpSACHGls4ZnVVVw5ZQQDcrOCjiMSGJWCCPDcup0camhmgU4dSZpTKYgQXTfh1IK+nH3q4KCjiARKpSBpb3tNHW+/t4/5xaO1upqkPZWCpL2ScISMkHHjLK2uJqJSkLTW1NLK42WVXDJhKEP75wYdRyRwKgVJa69trGZPXYMmmEU8KgVJayWlEYb2z+Hi8VpdTQRUCpLGdh88ymubqrlp1igyM/SjIAIqBUljj5VV0KrV1UQ+xLdSMLP7zazazMo72X6xmdWa2Wrv4wd+ZRE53rHV1c45bTBjC/oGHUckYfh5pPAAMOcEY/7gnJvuffzIxywiH/KXd/fx/t7DukS2yHF8KwXn3Apgn1/3L9ITy8MR+udmcsXkEUFHEUkoQc8pfNzM1pjZ82Y2qbNBZrbYzMJmFq6pqYlnPklBtUeaeG7dTq6frtXVRI4XZCmsBE5xzk0Dfg481dlA59xS51yxc664sFAvHZSeeXp1JQ3NWl1NpCOBlYJz7qBzrs67/RyQZWYFQeWR9FESjnDWiAFMLsoPOopIwgmsFMxsuHlXHzOzs70se4PKI+mhvLKW8sqD3Hy2jhJEOpLp1x2b2TLgYqDAzCqAO4EsAOfcEuAm4Mtm1gwcAW52zjm/8ohAdII5OzPEddN08TuRjvhWCs65hSfYfg9wj1/7Fzne0aYWnlpVyRWTh5Ofp9XVRDoS9KuPROLmhfJdHDzarAlmkS6oFCRtlJRGGDM4j3NOHRJ0FJGEpVKQtPD+3nr+vH0v84tHEQppdTWRzqgUJC0sD0cIGdw0S6eORLqiUpCU19zSymNlFVw8fijD87W6mkhXVAqS8t7YXMPugw2aYBaJgUpBUl5JaYSCfjlcMmFo0FFEEp5KQVJa9aGjvLqxmhtnFpGl1dVETkg/JZLSnlhZSXOrY75OHYnERKUgKcs5x/LSCLPHDuL0wn5BxxFJCioFSVnh9/ezfU+91mAW6QaVgqSsR96O0C8nk6umanU1kVipFCQlHToaXV3tmmkjycv27bqPIilHpSAp6Zk1OznS1KL3Joh0k0pBUlJJ6Q4mDO/PtFFaXU2kO1QKknI27jrImopa5hePxlvcT0RipFKQlFNSGiE7I8TcGVpdTaS7VAqSUhqaW3hyVSWXTxrGoL7ZQccRSToqBUkpL67fzYHDTZpgFjlJKgVJKSWlEYoG9uG80wuCjiKSlFQKkjIi+w7z5tY9zC8erdXVRE6SSkFSxqNlFZjBTcWjgo4ikrR8KwUzu9/Mqs2s/ATjZptZi5nd5FcWSX0trY5HwxEuHFdI0cA+QccRSVp+Hik8AMzpaoCZZQB3Ab/3MYekgT9sqWFn7VFNMIv0kG+l4JxbAew7wbCvAY8D1X7lkPSwPBxhcN9sLps4LOgoIkktsDkFMysC5gJLgsogqWFvXQMvvbObG2YUkZ2paTKRngjyJ+g/ge8551pONNDMFptZ2MzCNTU1cYgmyeTJVZU0tTidOhLpBUFeU7gYeMS7Nk0BcKWZNTvnnjp+oHNuKbAUoLi42MU1pSQ05xwlpRFmjBnIuGH9g44jkvQCKwXn3KnHbpvZA8DvOioEka6s3HGALdV13HXjlKCjiKQE30rBzJYBFwMFZlYB3AlkATjnNI8gvWJ5aYS87Ayumjoy6CgiKSHmUjCz84Fxzrn/b2aFQD/n3LudjXfOLYz1vp1zn411rMgxdQ3NPLO2iqunjqBfjlZXE+kNMU00m9mdwPeAv/M+lQX82q9QIrF4dm0VhxtbWDB7TNBRRFJGrK8+mgtcC9QDOOeqAM3qSaBKSiOcMbQfM8cMDDqKSMqItRQanXMOcABm1te/SCIntmX3IVbuOMACra4m0qtiLYXlZnYfMNDMvgS8DPzSv1giXSspjZCVYcydqdXVRHpTTLNzzrmfmNkngYPAeOAHzrmXfE0m0onG5laeWFXJZROHUdAvJ+g4IiklplLwThe96px7yczGA+PNLMs51+RvPJGPennDbvbVNzJf72AW6XWxnj5aAeR41yt6Gfgc0augisRdSWmEEfm5XDiuMOgoIikn1lIw59xh4Abg5865ucBZ/sUS6VjVgSOs2FLDvFmjyNDqaiK9LuZSMLOPA7cAz3qf07uFJO4eDVcAMK9Yp45E/BBrKXwd+D7whHNuvZmdCrzqXyyRj2ptdTxaFuG80wsYPTgv6DgiKSnWUjgMtAILzWwt8DTwCd9S+eDdPfV89TcrOdzYHHQUOUl/2raXiv1HNMEs4qNYTwE9DHwbKCdaDklnx77DPLduJw64Z+EMveEpCT1SuoOBeVlcfpZWVxPxS6xHCjXOuWecc+86594/9uFrsl520ZmFfOdTE3h27U6WvLE96DjSTfvrG3lx/W6un15EblZG0HFEUlasRwp3mtmvgFeAhmOfdM494Usqn9x+0WmUV9Xyb7/fyMQR/bl4/NCgI0mMnlpdSWNLq1ZXE/FZrEcKnwOmA3OAa7yPq/0K5Rcz4+6bpjJ+WH/+Ztkq3ttTH3QkicGx1dWmjcpn4ogBQccRSWmxlsI051yxc+4zzrnPeR+f9zWZT/KyM1m6qJhQyFj8UJj6Bk08J7q1FbVs3HVIE8wicRBrKbxlZinzZrUxQ/K4Z+FMtlbX8e1H1xC9AKwkqpJwhNysENdM0+pqIn6LtRTOB1ab2SYzW2tm67yXpiat88cV8HdXTOT58l3c+/q2oONIJw43NvP06iqumjKSAblZQccRSXmxTjTP8TVFQL54wamUV9Xykxc3MXFEfy6ZoJc6Jprn1u2irqFZE8wicRLTkUL7l6Em60tSO2Jm/OsNU5k4fABff2Q122vqgo4kx1leGuG0gr7MHjso6CgiaSHW00cpq092Bktvm0VWRojFD5Vx6KiuBp4ottXU8fZ7+5g/W6uricRL2pcCwKhBedzz6Rm8u6eeby1fQ2urJp4TwfJwhIyQcYNWVxOJG5WC59zTC7jjyom89M5ufv7q1qDjpL2mllYeL6vgkglDGdo/N+g4ImnDt1Iws/vNrNrMyjvZfp33SqbVZhY2s/P9yhKrz503lhtmFvHTlzfz0ju7g46T1l7dWM2eukZu1gSzSFz5eaTwAF2/aukVom+Kmw58HviVj1liYmb889wpTCnK55slq9larYnnoCwvjTC0fw4XnanV1UTiybdScM6tAPZ1sb3OffCusb5AQpzIz83K4L5Fs8jJDLH4wTAHNfEcd7tqj/LapmrmFY8iM0NnOEXiKdCfODOba2Ybia7mljCXzRg5sA/33jKTHfsO881HVmviOc4eX1lBq4P5Wl1NJO4CLQXn3JPOuQnA9cCPOxtnZou9eYdwTU1NXLJ97LQh/OCas3hlYzX/+fLmuOxToqurLQ9HOOe0wZwypG/QcUTSTkIcm3unmk43s4JOti/1LshXXFgYv3PMi845hXmzRvGzV7fyQvmuuO03nb317l7e33uYm2ePCTqKSFoKrBTM7Azz3pFkZjOBbGBvUHk6Ymb8+PrJTBs9kL9dvpotuw8FHSnlLS+N0D83kzmThwcdRSQt+fmS1GXAn4HxZlZhZl8ws9vN7HZvyI1AuZmtBn4BLHAJeLnS3KwM7rt1Fn2yM/nSg2Fqj2ji2S+1h5t4vnyXVlcTCVCsF8TrNufcwhNsvwu4y6/996bh+bksuXUmC3/5Fl9/ZBX/85nZZIR02YXe9ts1lTQ0a3U1kSAlxJxCMigeO5g7r5nE65tq+PcXNwUdJyWVlEaYNHIAk4vyg44ikrZUCt1wy8fGsPDs0dz7+jaeXbsz6DgppbyylvVVB3WUIBIwlUI3mBk/vHYSM8cM5NuPrmHjroNBR0oZJaURcjJDXDdNF78TCZJKoZtyMjNYcuss+udmsvjBMg4cbgw6UtI72tTCU6sruWLycPLztLqaSJBUCidh6IBcliyaxa7ao3xt2Spa9I7nHnmhfBeHjjYzX6eORAKnUjhJM8cM4kfXTeIPW/bwb7/fGHScpPZI6Q5OGZLHOacOCTqKSNpTKfTAzWeP4dZzxnDfG9t5ek1V0HGS0nt76nlr+z7mF48mpJf5igROpdBDP7h6ErPHDuK7j63hnSpNPHfX8nCEkMFNs0YFHUVEUCn0WHZmiF/cMpOBfbJZ/FCYffWaeI5Vc0srj5VV8InxQxk2QKuriSQClUIvGNo/OvFcfaiBr/5mJc0trUFHSgpvbK6h+lCDJphFEohKoZdMHz2Qf7p+Mn/atpd/fV4Tz7F4pDRCQb8cLpkwNOgoIuJRKfSi+cWj+ey5Y/nVm+/y5KqKoOMktOpDR3l1YzU3zioiS6uriSQM/TT2sjuumsjHTh3M9x9fR3llbdBxEtYTKytpaXVaXU0kwagUellWRnTieUjfbP7qoTL21jUEHSnhOOdYXhrh7LGDOb2wX9BxRKQdlYIPCvrlcN+iYvbUNfCV36ykSRPPH1L63n6276nXBLNIAlIp+GTKqHz+5YYpvLV9H//v2Q1Bx0koJaUR+uVkcuUUra4mkmh8W2RH4IaZoyivPMj9f3yXyUX5eoMWcPBoE8+uq+KGmaPIy9Z/P5FEoyMFn/39lRM49/Qh/P2T61gTORB0nMA9s6aKo02tLNAEs0hCUin4LDMjxD2fnklhvxz+6qEyag6l98Tz8tIIE4b3Z+oora4mkohUCnEwuG829y2axYEjjfz1w2U0NqfnxPOGnQdZU1HLgtmjMdPF70QSkUohTiYX5XPXjVMpfW8/P/7dO0HHCURJaYTsjBDXT9fqaiKJSjN9cXTd9CLWVx1k6YrtTC4awILZY4KOFDfHVle7fNIwBvXNDjqOiHRCRwpx9t1PjeeCcQX8w1PrWbljf9Bx4ubFd3Zz4HATN6dREYokI99KwczuN7NqMyvvZPstZrbW+/iTmU3zK0siycwI8fOFMxiWn8OXf11G9cGjQUeKi+WlEUYN6sO5p2t1NZFE5ueRwgPAnC62vwtc5JybCvwYWOpjloQyMC+bpYuKOXikmS8/vDLlJ54j+w7z5tY9Wl1NJAn4VgrOuRXAvi62/8k5d+z8yVtAWr2za+KIAdw9bypl7+/nh8+sDzqOrx4NRzCtriaSFBJlovkLwPNBh4i3q6eOZH3VQf779W1MGjmAWz52StCRel1Lq+PRsgouHFfIyIF9go4jIicQ+ESzmX2CaCl8r4sxi80sbGbhmpqa+IWLg29fPp6Lzizkh0+vJ/xepwdWSWvFlhp21h7lZl38TiQpBFoKZjYV+BVwnXNub2fjnHNLnXPFzrniwsLC+AWMg4yQ8bObZzByYB++/PBKdtWm1sTz8tIIQ/pmc+nEYUFHEZEYBFYKZjYGeAJY5JzbHFSORJCfl8UvbyumvqGZ239dRkNzS9CResWeugZe3rCbuTOKyM4M/KBURGLg50tSlwF/BsabWYWZfcHMbjez270hPwCGAPea2WozC/uVJRmcOaw//zF/GqsjB/iHp8pxzgUdqceeXFlJU4tjgU4diSQN3yaanXMLT7D9i8AX/dp/MpozeQRfu+QMfv7qVqYU5bPo42ODjnTSnHOUhCPMHDOQccP6Bx1HRGKkY/oE883LzuSSCUP5x2fe4e13k3fieeWOA2ytrtNRgkiSUSkkmFDI+OmC6YwZnMdfP1xG1YEjQUc6KSWlO+ibncHVU0cGHUVEukGlkIDy+2Sx9LZZHG1q5fZfl3G0Kbkmnusamvnd2p1cPXUkfXMS5a0wIhILlUKCOmNodOJ5bUUtdzyZXBPPz66t4nBjC/N16kgk6agUEtjlk4bz9UvH8fjKCv73T+8FHSdmj5RGGDe0HzPHDAw6ioh0k0ohwX390nFcNnEYP352A3/e1un7+xLGlt2HWLXjgFZXE0lSKoUEF514nsbYIXl85Tcrqdh/OOhIXSopjZCVYcydodXVRJKRSiEJ9M/NYultxTQ1t/JXD5VxpDExJ54bm1t5YlUlnzxrGEP65QQdR0ROgkohSZxe2I//Wjidd3Ye5O+eWJuQE88vb9jNvvpG5hdrglkkWakUksglE4bxrcvO5KnVVfzPm+8GHecjSkojjMzP5YJxqXXRQpF0olJIMl/5xBnMmTScf35uA3/cuifoOG0qDxxhxZYabioeTYZWVxNJWiqFJBMKGT+ZP43TC/vx1d+sJLIvMSaeHwtXADBPq6uJJDWVQhLql5PJL28rpqXVsTgBJp5bWx3LwxHOP6OA0YPzAs0iIj2jUkhSYwv68rOFM9i46yDffTzYiec/bttD5YEjmmAWSQEqhSR28fihfOdT43lmTRVLV2wPLEdJaYSBeVlcPkmrq4kkO5VCkvvyRadz1ZQR3PXCRlZsjv/61fvrG3lxfXR1tZzMjLjvX0R6l0ohyZkZd8+bypnD+vO1Zat4f299XPf/5KpKGltatW6CSIpQKaSAvOxMli4qBmDxg2XUNzTHZb/OOUpKI0wblc+E4QPisk8R8ZdKIUWMGZLHPZ+ewZbqQ3znsTVxmXheU1HLpt2HWDB7jO/7EpH4UCmkkAvGFfL9Kybw3Lpd3Pv6Nt/3V1IaoU9WBtdMG+H7vkQkPlQKKeZLF5zGtdNG8pMXN/Hapmrf9nO4sZln1lRx5ZQR9M/N8m0/IhJfKoUUY2bcdeNUJg4fwN8sW8W7e/yZeH527U7qGpq5+WxNMIukEpVCCuqTncF9i2aRGTIWPximzoeJ5+XhCKcV9qX4lEG9ft8iEhzfSsHM7jezajMr72T7BDP7s5k1mNm3/cqRrkYPzuOeT89k+556vlWymtbW3pt43lZTR+l7+1lQrNXVRFKNn0cKDwBzuti+D/gb4Cc+Zkhr551RwN9fOZEX39nNPa9t7bX7XV4aITNk3DBTF78TSTW+lYJzbgXRJ/7Otlc750qBJr8yCHz+vLHMnVHET1/ezMvv7O7x/TW1tPL4ygoumTCUwv5aXU0k1STFnIKZLTazsJmFa2rifymHZGZm/MsNU5g0cgDfLFnNtpq6Ht3fqxur2VPXqAlmkRSVFKXgnFvqnCt2zhUXFmpVr+7KzcrgvkXFZGeG+NKDYQ4ePfmDs5LSCMMG5HChVlcTSUlJUQrSc0UD+/CLW2by/t7DJz3xvKv2KK9vquamWaPIzNB/HZFUpJ/sNHLOaUP4h6sm8vKGav7rlS3d/vrHyiK0OrRugkgKy/Trjs1sGXAxUGBmFcCdQBaAc26JmQ0HwsAAoNXMvgGc5Zw76Fcmgc+cO5byqoP81ytbOGvkAD41aXhMXxddXa2Cj582hFOG9PU5pYgExbdScM4tPMH2XYBe0xhnZsY/XT+ZLbsP8a2S1Tz1lfMYN6z/Cb/urXf3smPfYb71yTPjkFJEgqLTR2koNyuDJYtm0Sc7g8UPlVF75MQTzyWlEQbkZjJncmxHFiKSnFQKaWpEfh/++9ZZRPYd5huPrKKli4nn2sNNPF++i+tnFJGbpdXVRFKZSiGNzR47mDuvncRrm2r46UubOx332zWVNDa3aoJZJA34NqcgyeHWj41hfWUt97y2lUkjB3DFlI+ujfDI2xEmFw1gclF+AAlFJJ50pJDmzIx/vG4SM8YM5G8fXcOmXYc+tL28spZ3dh5kgY4SRNKCSkHIycxgya2z6JuTyeKHwtQe/mDiuaQ0Qk5miGunFwWYUETiRaUgAAwbkMuSW2dSdeAIX/Mmno82tfDU6kqunDKC/D5aXU0kHagUpM2sUwbzo+sms2JzDXf/fhPPl+/k0NFmTTCLpBFNNMuHLDx7DOsqa1nyxjaGD8hl7JA8zjltcNCxRCROdKQgH/HDayZRfMogdh08yjytriaSVnSkIB+RnRni3ltn8sAf3+PWc04JOo6IxJFKQTo0tH8u350zIegYIhJnOn0kIiJtVAoiItJGpSAiIm1UCiIi0kalICIibVQKIiLSRqUgIiJtVAoiItLGnOt8GcZEZGY1wPsn+eUFwJ5ejNNbEjUXJG425eoe5eqeVMx1inOu8ESDkq4UesLMws654qBzHC9Rc0HiZlOu7lGu7knnXDp9JCIibVQKIiLSJt1KYWnQATqRqLkgcbMpV/coV/ekba60mlMQEZGupduRgoiIdCGlS8HM7jazjWa21syeNLOBnYybY2abzGyrmX0/Drnmmdl6M2s1s05fSWBm75nZOjNbbWbhBMoV18fL2+dgM3vJzLZ4fw7qZFyL93itNrOnfcrS5fdvZjlmVuJt/4uZjfUjx0nk+qyZ1bR7fL4Yp1z3m1m1mZV3st3M7Gde7rVmNjNBcl1sZrXtHq8fxCnXaDN7zcw2eD+PX+9gjH+PmXMuZT+Ay4FM7/ZdwF0djMkAtgGnAdnAGuAsn3NNBMYDrwPFXYx7DyiI4+N1wlxBPF7efv8N+L53+/sd/Vt62+p8znHC7x/4a2CJd/tmoCQOj08suT4L3BOv/0/t9nshMBMo72T7lcDzgAHnAH9JkFwXA78L4PEaAcz0bvcHNnfwb+nbY5bSRwrOuRedc83eX98CRnUw7Gxgq3Nuu3OuEXgEuM7nXBucc5v83MfJiDFX3B8vz3XA/3q3/xe4Pg777Egs33/7rI8Bl5r/C10H9e9yQs65FcC+LoZcBzzoot4CBprZiATIFQjn3E7n3Erv9iFgA1B03DDfHrOULoXjfJ5osx6vCIi0+3sFH/0HCIoDXjSzMjNbHHQYT1CP1zDn3E6I/tAAQzsZl2tmYTN7y8z8KI5Yvv+2Md4vJbXAEB+XtmT/AAAFVklEQVSydDcXwI3e6YbHzGy0z5lilcg/gx83szVm9ryZTYr3zr1TjzOAvxy3ybfHLOnXaDazl4HhHWy6wzn3W2/MHUAz8HBHd9HB53r8kqxYcsXgPOdclZkNBV4ys43ebzdB5vLl8YKus3XjbsZ4j9lpwKtmts45t6038nli+f59e4y6EMs+nwGWOecazOx2okczl/icKxZBPF6xWEn00hB1ZnYl8BQwLl47N7N+wOPAN5xzB4/f3MGX9MpjlvSl4Jy7rKvtZvYZ4GrgUuedjDtOBdD+N6ZRQJXfuWK8jyrvz2oze5LoKYIelUIv5PLl8YKus5nZbjMb4Zzb6R0mV3dyH8ces+1m9jrR37J6sxRi+f6Pjakws0wgH/9PU5wwl3Nub7u//pLoPFsi8O3/VE+0fyJ2zj1nZveaWYFzzvdrIplZFtFCeNg590QHQ3x7zFL69JGZzQG+B1zrnDvcybBSYJyZnWpm2UQnBn151Up3mFlfM+t/7DbRSfMOXyURZ0E9Xk8Dn/Fufwb4yFGNmQ0ysxzvdgFwHvBOL+eI5ftvn/Um4NVOfiGJa67jzjlfS/RcdSJ4GrjNe0XNOUDtsVOFQTKz4cfmgszsbKLPl3u7/qpe2a8B/wNscM79RyfD/HvM4j2zHs8PYCvR826rvY9jrwgZCTzXbtyVRGf4txE9jeJ3rrlEm74B2A38/vhcRF9Fssb7WJ8ouYJ4vLx9DgFeAbZ4fw72Pl8M/Mq7fS6wznvM1gFf8CnLR75/4EdEf/kAyAUe9f7/vQ2cFqfH6ES5/sX7v7QGeA2YEKdcy4CdQJP3/+sLwO3A7d52A37h5V5HF6/Ii3Our7Z7vN4Czo1TrvOJngpa2+6568p4PWZ6R7OIiLRJ6dNHIiLSPSoFERFpo1IQEZE2KgUREWmjUhARkTYqBUkbZlbXw69/zHundFdjXrcurjAb65jjxhea2QuxjhfpCZWCSAy8695kOOe2x3vfzrkaYKeZnRfvfUv6USlI2vHeBXq3mZVbdL2KBd7nQ96lDNab2e/M7Dkzu8n7slto9y5qM/tv78J7683sHzvZT52Z/buZrTSzV8yssN3meWb2tpltNrMLvPFjzewP3viVZnZuu/FPeRlEfKVSkHR0AzAdmAZcBtztXQLiBmAsMAX4IvDxdl9zHlDW7u93OOeKganARWY2tYP99AVWOudmAm8Ad7bblumcOxv4RrvPVwOf9MYvAH7WbnwYuKD736pI9yT9BfFETsL5RK8W2gLsNrM3gNne5x91zrUCu8zstXZfMwKoaff3+d7lzDO9bWcRvSxBe61AiXf710D7C5sdu11GtIgAsoB7zGw60AKc2W58NdHLjYj4SqUg6aizBW+6WgjnCNFrGmFmpwLfBmY75/ab2QPHtp1A+2vKNHh/tvDBz+E3iV5zahrRo/ij7cbnehlEfKXTR5KOVgALzCzDO89/IdEL171JdBGakJkNI7oc4zEbgDO82wOAeqDWG3dFJ/sJEb1KKsCnvfvvSj6w0ztSWUR0ic1jziQxrpIrKU5HCpKOniQ6X7CG6G/v33XO7TKzx4FLiT75bia62lWt9zXPEi2Jl51za8xsFdEraG4H/tjJfuqBSWZW5t3PghPkuhd43MzmEb2KaX27bZ/wMoj4SldJFWnHzPq56EpbQ4gePZznFUYfok/U53lzEbHcV51zrl8v5VoBXOec298b9yfSGR0piHzY78xsIJAN/Ng5twvAOXfEzO4kug7ujngG8k5x/YcKQeJBRwoiItJGE80iItJGpSAiIm1UCiIi0kalICIibVQKIiLSRqUgIiJt/g+3apY7akfpKwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3ce4b88390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('alpha is:', 0.1)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef_lasso</th>\n",
       "      <th>coef_lr</th>\n",
       "      <th>coef_ridge</th>\n",
       "      <th>columns</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.235076</td>\n",
       "      <td>[0.672673637912808]</td>\n",
       "      <td>[0.3483407201348778]</td>\n",
       "      <td>temp</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>[0.5814430328088782]</td>\n",
       "      <td>[0.07049018062026347]</td>\n",
       "      <td>mnth</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.223094</td>\n",
       "      <td>[0.2928577641907567]</td>\n",
       "      <td>[0.27883449919778436]</td>\n",
       "      <td>season</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>[0.03660845226806475]</td>\n",
       "      <td>[0.03574652365761666]</td>\n",
       "      <td>weekday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>[0.005392332572906798]</td>\n",
       "      <td>[0.004872538193709086]</td>\n",
       "      <td>workingday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>[6.106226635438361e-16]</td>\n",
       "      <td>[0.0]</td>\n",
       "      <td>yr</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-0.000000</td>\n",
       "      <td>[-0.036070039662185804]</td>\n",
       "      <td>[-0.03427217730207488]</td>\n",
       "      <td>holiday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.335821</td>\n",
       "      <td>[-0.03642170421782992]</td>\n",
       "      <td>[0.2842833498135113]</td>\n",
       "      <td>atemp</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>-0.000000</td>\n",
       "      <td>[-0.05659929339172948]</td>\n",
       "      <td>[-0.06537306242474539]</td>\n",
       "      <td>hum</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>-0.026223</td>\n",
       "      <td>[-0.11474730632971941]</td>\n",
       "      <td>[-0.10957710086690556]</td>\n",
       "      <td>windspeed</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>-0.171476</td>\n",
       "      <td>[-0.23037733832290053]</td>\n",
       "      <td>[-0.21658061667954276]</td>\n",
       "      <td>weathersit</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>[-0.5785162630911328]</td>\n",
       "      <td>[-0.05613399264135427]</td>\n",
       "      <td>instant</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    coef_lasso                  coef_lr              coef_ridge     columns\n",
       "8     0.235076      [0.672673637912808]    [0.3483407201348778]        temp\n",
       "3     0.000000     [0.5814430328088782]   [0.07049018062026347]        mnth\n",
       "1     0.223094     [0.2928577641907567]   [0.27883449919778436]      season\n",
       "5     0.000000    [0.03660845226806475]   [0.03574652365761666]     weekday\n",
       "6     0.000000   [0.005392332572906798]  [0.004872538193709086]  workingday\n",
       "2     0.000000  [6.106226635438361e-16]                   [0.0]          yr\n",
       "4    -0.000000  [-0.036070039662185804]  [-0.03427217730207488]     holiday\n",
       "9     0.335821   [-0.03642170421782992]    [0.2842833498135113]       atemp\n",
       "10   -0.000000   [-0.05659929339172948]  [-0.06537306242474539]         hum\n",
       "11   -0.026223   [-0.11474730632971941]  [-0.10957710086690556]   windspeed\n",
       "7    -0.171476   [-0.23037733832290053]  [-0.21658061667954276]  weathersit\n",
       "0     0.000000    [-0.5785162630911328]  [-0.05613399264135427]     instant"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mses = np.mean(lasso.mse_path_, axis = 1)\n",
    "plt.plot(np.log10(lasso.alphas_), mses) \n",
    "#plt.plot(np.log10(lasso.alphas_)*np.ones(3), [0.3, 0.4, 1.0])\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show()    \n",
    "            \n",
    "print ('alpha is:', lasso.alpha_)\n",
    "\n",
    "# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性\n",
    "fs = pd.DataFrame({\"columns\":list(columns), \"coef_lr\":list((lr.coef_.T)), \"coef_ridge\":list((ridge.coef_.T)), \"coef_lasso\":list((lasso.coef_.T))})\n",
    "fs.sort_values(by=['coef_lr'],ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
